System identification method

ABSTRACT

Quick real-time identification and estimation of a time-non-varying or time-varying system. A new H ∞  evaluation criterion is determined, a fast algorithm for a modified H ∞  filter based on the criterion is developed, and a quick time-varying system identifying method according to the fast H ∞  filtering algorithm is provided. By the fast H ∞  filtering algorithm, a time-varying system sharply varying can be traced with an amount of calculation O(N) per unit time step. The algorithm completely agrees with a fast Kalman filtering algorithm at the extreme of the upper limit value. If the estimate of impulse response is determined, a pseudo-echo is sequentially determined from the estimate and subtracted from the real echo to cancel the echo. Thus an echo chancellor is realized.

CROSS REFERENCES TO RELATED APPLICATIONS

This application is based upon priority International ApplicationPCT/JP01/09082 filed Oct. 16, 2001, International Publication No. WO02/35727 A1 published May 2, 2002, which is based upon JapaneseApplication 2000-323958 filed Oct. 24, 2000.

BACKGROUND OF THE INVENTION

The present invention relates to system identification methods, and moreparticularly, to a time-varying system identification method at a highspeed in real time by using a fast algorithm for modified H_(∞) filtersdeveloped based on a new H_(∞) evaluation criterion.

In general, system identification is to estimate a mathematical model (atransfer function, an impulse response, or the like) of a systeminput-and-output relationship according to the input-and-output data.Typical applications thereof include echo cancellers in internationalcommunications, automatic equalizers in data communications, echocancellers in acoustic systems, sound-field reproduction, and activenoise control in vehicles. Details are written in “Digital SignalProcessing Handbook”, the Institute of Electronics, Information andCommunication Engineers, 1993, or the like.

(Basic Principle)

FIG. 14 shows a configuration for system identification. This systemincludes an unknown system 1 and an adaptive filter 2. The adaptivefilter 2 has an FIR digital filter 3 and an adaptive algorithm 4.

One case which uses an output error method to identify the unknownsystem 1 will be described below. Here, u_(k) indicates an input to theunknown system 1, d_(k) indicates the output of the system, which is asignal to be obtained, and d^_(k) indicates the output of the filter. (Amark “^” means an estimated value and should be placed above characters,but it is placed at the upper right of the characters for inputconvenience. This notation may be used through the presentspecification.)

Since an impulse response is generally used as a parameter of an unknownsystem, the adaptive filter adjusts a coefficient of the FIR digitalfilter 3 by the adaptive algorithm such that an evaluation errore_(k)=d_(k)−d^_(k) in the figure is minimized.

FIG. 15 shows the structure of an impulse-response adjustment mechanism.

Here, as an example of adaptive algorithm, the following LMS algorithm(least mean square algorithm) is widely used because of itscomputational simplicity.

[LMS Algorithm]ĥ _(k+1) =ĥ _(k) +μu _(k)(y _(k) −u _(k) ^(T) ĥ _(k))  (1)where,ĥ _(k) =[ĥ ₀ [k], . . . ,ĥ _(N−1) [k]] ^(T) , u _(k) =[u _(k) , . . . ,u _(k−N+1)]^(T), μ>0  (2)Generally, Kalman filters, which converges quickly, are suitable foridentifying a time-varying system.[Kalman Filter]{circumflex over (x)} _(k|k) ={circumflex over (x)} _(k|k−1) +K _(k)(y_(k) −H _(k) {circumflex over (x)} _(k|k−1)){circumflex over (x)}_(k+1|k)={circumflex over (x)}_(k|k)  (3)K _(k) ={circumflex over (P)} _(k|k−1) H _(k) ^(T)(1+H _(k) {circumflexover (P)} _(k|k−1) H _(k) ^(T))⁻¹  (4){circumflex over (P)} _(k|k) ={circumflex over (P)} _(k|k−1) −K _(k) H_(k) {circumflex over (P)} _(k|k−1)

$\begin{matrix}{{\hat{P}}_{{k + 1}|k} = {{\hat{P}}_{k|k} + {\frac{\sigma_{\omega}^{2}}{\sigma_{\upsilon}^{2}}I}}} & (5)\end{matrix}$where,{circumflex over (x)} _(k|k) =[ĥ ₀ [k], . . . , ĥ _(N−1) [k]] ^(T) , H_(k) =[u _(k−1) , . . . , u _(k−N]){circumflex over (x)} _(0|−1)=0, {circumflex over (P)} _(0|−1)=ε₀ I, ε₀>0  (6)Here, the impulse response {h_(i)} of unknown system is obtained as astate estimate x^_(k|k), and an input {u_(k)} to the system is used asan element of an observation matrix H_(k).

Fast Kalman filtering algorithm is also known, which reduces the amountof calculation per unit time step to the number of times calculationsare performed proportional to N, that is, O(N) by applying the shiftproperty (H_(k+1)(i+1)=H_(k)(i)) of the observation matrix H_(k) to aKalman filter obtained when σ² _(w)=0. Details are written in “DigitalSignal Processing Handbook”, the Institute of Electronics, Informationand Communication Engineers, 1993, or the like.

(Applications to Echo Cancellers)

Four-wire circuits are used for long distance telephone lines such asfor international calls for a reason of signal amplification and others.On the other hand, two-wire circuits are used for subscriber linesbecause they are relatively short.

FIG. 16 is an explanation view of a communication system and an echo.Impedance matching is performed by disposing hybrid transformers atconnection points between two-wire circuits and a four-wire circuit, asshown in the figure. If the impedance matching is perfect, a signal(voice) from a speaker B reaches only a speaker A. However, it isdifficult to make the matching perfect in general. A phenomenon occursin which a part of a received signal leaks to the four-wire circuit, isamplified, and returns to the receiver (speaker A). This is an echo. Asthe transmission distance gets longer (delay time gets longer), theeffect of the echo gets larger, and the quality of telephone speechsignificantly deteriorates (in pulse transmission, an echo influencessignificantly even in a short distance, and the quality of telephonespeech deteriorates). FIG. 17 shows a basic principle of an echocanceller.

As shown in the figure, the echo canceller is introduced to successivelyestimate the impulse response of an echo path by using a received signaland an echo directly observable, and to subtract a quasi echo obtainedby using the estimate from an actual echo to cancel and remove the echo.

The impulse response of the echo path is estimated such that the meansquare error of a residual echo e_(k) is minimized. Elements whichdisturb the estimation of the echo path are line noise and a signal(voice) from the speaker A. When two speakers start talking at the sametime (double talking), the estimation of the impulse response isgenerally halted. Since the impulse response of the hybrid transformershas a length of approximately 50 [ms], if the sampling period is set to125 [μs], the order of the impulse response of the echo path actuallybecomes approximately 400.

BRIEF SUMMARY OF THE INVENTION

In conventional arts, the LMS algorithm (least-mean-square algorithm)has been widely used as adaptive algorithm due to its simplicity, but itis impossible to closely track a time-varying system which variessuddenly due to its very slow convergence.

A Kalman filter, which has an excellent tracking performance, requiresthe amount of calculation proportional to O(N²) or O(N³). Since theamount of calculation rapidly increases with the tap number N, it isdifficult to process actual issues requiring a high tap number N in realtime. As a countermeasure therefor, a fast Kalman filter with thecomputational complexity of O(N) per unit time step for a tap number Nhas been proposed, but it is impossible for the filter to identify atime-varying system because of its stationary characteristic (incapableof taking system noise into account).

In view of the above-described points, it is an object of the presentinvention to implement a fast real-time identification of time-varyingand time-invariant systems by using a fast algorithm for modified H_(∞)filters developed based on a new H_(∞) evaluation criterion. It isanother object of the present invention to include, as a particular caseof the present algorithm, a fast Kalman filtering algorithm, and todetermine theoretically the covariance of system noise which is dominantin the tracking performance of time-varying systems. It is still anotherobject of the present invention to provide a fast time-varying systemidentification method which is substantially effective even when aninput signal is discontinuously varied, such as in an echo canceller fora time-varying system which varies extremely as sudden line switching.It is a further object of the present invention to provide a systemidentification method which is applicable to echo cancellers incommunication systems and acoustic systems, sound-field reproduction,and noise control.

In order to solve the problems described above, a new H_(∞) evaluationcriterion is introduced, a fast algorithm for the modified H_(∞) filtersis developed according to the criterion, and a fast time-varying systemidentification method based on the fast algorithm is proposed in thepresent invention. The fast algorithm according to the present inventionis capable of tracking a time-varying system which varies rapidly, withthe complexity of O(N) per unit time step. Further, it has a convenientproperty that it completely matches the fast Kalman filtering algorithmat a limit of γ_(f)=∞.

A more detailed explanation of the invention is provided in thefollowing description and appended claims take in conjunction with theaccompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flowchart of a fast algorithm.

FIG. 2 is an explanation view (1) of the amount of calculation in eachpart of a modified H_(∞) filtering algorithm.

FIG. 3 is an explanation view (2) of the amount of calculation in eachpart of the modified H_(∞) filtering algorithm.

FIG. 4 is an explanation view of the amount of calculation in theRiccati recursion.

FIG. 5 is an explanation view (1) of the amount of calculation in a fastH_(∞) filtering algorithm.

FIG. 6 is an explanation view (2) of the amount of calculation in thefast H_(∞) filtering algorithm.

FIG. 7 is a view showing values of an impulse response {h_(i)} to beestimated.

FIG. 8 is a comparative explanation view of estimated results of impulseresponses obtained by the modified H_(∞) filter and the fast H_(∞)filter.

FIG. 9 is a view of measurement results of computation time.

FIG. 10 is a view (1) of simulation results of each algorithm.

FIG. 11 is a view (2) of simulation results of each algorithm.

FIG. 12 is a view showing the relationship between γ_(f) and ρ.

FIG. 13 is a view showing the relationship among tap numbers andcomputation times of impulse responses obtained by the fast H_(∞)filter, fast Kalman filter, and LMS algorithm.

FIG. 14 is a view showing a configuration for system identification.

FIG. 15 is a view showing the configuration of an impulse-responseadjustment mechanism.

FIG. 16 is an explanation view of a communication system and an echo.

FIG. 17 is a view showing the principle of an echo canceller.

DETAILED DESCRIPTION OF THE INVENTION

A detailed description of the preferred embodiments and best modes ofpracticing the present invention are described hereinafter. Details areshown, for example, in “Derivation of A Fast Algorithm of Modified H_(∞)Filters”, K. Nishiyama, IEEE international Conference on IndustrialElectronics, Control and Instrumentation, RBC-II, pp. 462–467, October,2000.

1. Description of Symbols

First, main symbols used in the embodiments of the present invention andwhether they are known or unknown will be described.

-   x_(k): State vector or just state, unknown and to be estimated.-   x₀: Initial state, unknown.-   w_(k): System noise, unknown.-   v_(k): Observation noise, unknown.-   y_(k): Observation signal, known and input to a filter.-   z_(k): Output signal, unknown.-   H_(k): Observation matrix, known.-   L_(k): Output matrix, known.-   x^_(k|k): State value of the state x_(k) at time k, estimated by    using observation signals y₀−y_(k). Given by a filter equation.-   x^_(0|0): Initial estimate of a state, essentially unknown but set    to 0 for convenience.-   K_(s,k+1): Filter gain, obtained by matrix P^_(k|k−1.)-   Σ_(wk): Corresponds to the covariance matrix of the system noise,    known in theory but unknown in practice.-   P^_(k|k−1): Corresponds to the covariance matrix of the error of    X^_(k|k−1), given by a Riccati equation.-   P^_(1|0): Corresponds to the covariance matrix of an error in the    initial state, essentially unknown but set to ε₀I for convenience.-   σ² _(v): Variance of the observation noise, treated as known in    theory but unknown in practice.-   σ² _(w): Variance of the system noise, treated as known in theory    but unknown in practice.

A mark “^” placed above a symbol indicates an estimated value, a mark“U” indicates that the matrix is extended by one row, and a mark “˜” isadded for convenience. These marks are placed at the upper right ofcharacters for input convenience, but, as shown in expressions, they areidentical with those placed above characters. “L”, “H”, “P”, and “K”indicate matrixes. Some of them are written in boldface as inexpressions, but they are usually written in lightface for convenience.

2. Modified H_(∞) Filter

Next, a state-space model as in the following equations (7) to (9) isdiscussed.x _(k+1) =x _(k) +w _(k) , w _(k) , x _(k) εR ^(N)  (7)y _(k) =H _(k) x _(k) +v _(k) , y _(k) , v _(k) εR  (8)z_(k)=H_(k)x_(k), z_(k)εR, H_(k)εR^(1×N)  (9)where, L_(k)=H_(k)(H_(k)=[u_(k)u_(k−2) . . . u_(k−N+1)]) assuming anecho canceller or the like. For such a state-space model, an H_(∞)evaluation criterion (γ_(f) is newly placed in the left-hand side) suchas that shown by expression (10) is proposed.

$\begin{matrix}{{\sup\limits_{x_{0},{\{\omega_{i}\}},{\{\upsilon_{i}\}}}\frac{\sum\limits_{i = 0}^{k}{{e_{f,i}}^{2}/\rho}}{{{x_{0} - {\overset{\Cup}{x}}_{0|{- 1}}}}_{\sum_{0}^{- 1}}^{2} + {\sum\limits_{i = 0}^{k}{\omega_{i}}_{\sum_{\omega_{i}}^{- 1}}^{2}} + {\sum\limits_{i = 0}^{k}{{\upsilon_{i}}^{2}/\rho}}}} < \gamma_{f}^{2}} & (10)\end{matrix}$

When it is assumed that ρ or Σ_(wk) does not depend on γ_(f), a modifiedH_(∞) filter of level γ_(f) satisfying the evaluation criterion can begiven by the following equations (11) to (14) by applying a standardH_(∞) estimation scheme known in the system identification field. Thisscheme is shown, for example, in “Linear Estimation in Krein Spaces—PartI: Theory,” B. Hassibi, A. H. Sayed, and T. Kailath, IEEE Trans.Automatic Control, 41, 1, pp. 18–33, 1996., and “Linear Estimation inKrein Spaces—Part II: Applications,” B. Hassibi, A. H. Sayed, and T.Kailath, IEEE Trans. Automatic Control, 41, 1, pp. 34–49, 1996.{hacek over (z)}_(k|k)=H_(k){circumflex over (x)}_(k|k)  (11){circumflex over (x)} _(k+1|k+1) ={circumflex over (x)} _(k|k) +K_(s,k+1)(y _(k+1) −H _(k+1) {circumflex over (x)} _(k|k)) Filterequation  (12)K _(s,k+1) ={circumflex over (P)} _(k+1|k) H _(k+1) ^(T)(H _(k+1){circumflex over (P)} _(k+1|k) H _(k+1) ^(T)+ρ)⁻¹ Filter gain  (13)

                            Riccati  equation  (14)${\hat{P}}_{k + {1{k}}} = {{\hat{P}}_{k{{k - 1}}} - {{{\hat{P}}_{k{{k - 1}}}\left\lbrack {H_{k}^{T}H_{k}^{T}} \right\rbrack}{R_{e,k}^{- 1}\begin{bmatrix}H_{k} \\H_{k}\end{bmatrix}}{\hat{P}}_{k{{k - 1}}}} + \Sigma_{wk}}$where,

$\begin{matrix}{{e_{j,i} = {{\overset{˘}{z}}_{i{i}} - {H_{i}x_{i}}}}{R_{e,k} = {R_{k} + {\begin{bmatrix}H_{k} \\H_{k}\end{bmatrix}{{\hat{P}}_{k{{k - 1}}}\left\lbrack {H_{k}^{T}H_{k}^{T}} \right\rbrack}}}}{{R_{k} = \begin{bmatrix}\rho & 0 \\0 & {- {\rho\gamma}_{f}^{2}}\end{bmatrix}},{\Sigma_{wk} = {\gamma_{f}^{- 2}{\hat{P}}_{k + {1{k}}}}}}{{{{\hat{P}}_{k{{k - 1}}}^{- 1} + {H_{k}^{T}H_{k}}} > 0},{{\hat{P}}_{1{0}} = {ɛ_{0}I}},{ɛ_{0} > 0}}{{{0 < \rho} = {{1 - \gamma_{f}^{- 2}} \leq 1}},{\gamma_{f} > 1}}} & (15)\end{matrix}$

Since a weight ρ in the evaluation criterion depends on an upper limitγ_(f) determined in advance, the above algorithm is essentiallydifferent from that applied to normal H_(∞) filters. The presentalgorithm controls a maximum energy gain from disturbances (having theinitial state x₀, the system noise {w_(i)}, and the observation noise{υ_(i)}) weighted by ρ to a filter error {e_(f,i)} so as to be smallerthan γ_(f) ². Therefore, the present algorithm is a robust filteringalgorithm against the disturbances. This property is reflected by thetracking characteristic of a time-varying system. When γ_(f)→∞ issatisfied, ρ=1 and Σ_(wk)=0. In this time, the modified H_(∞) filterbecomes a standard Kalman filter.

The main load for calculating the modified H_(∞) filter rises during theupdate of P^_(k+1|k)εR^(N×N), which requires the amount of calculationin proportion to N² or N³. That is, an arithmetic operation of O(N²) perunit time step is required. Here, a tap number N matches the dimensionof the state vector x_(k). Therefore, as the dimension of x_(k)increases, the computation time required to perform the modified H_(∞)filter increases rapidly. In order to overcome the drawback, theintroduction of a fast algorithm of the modified H_(∞) filter is needed.

3. Fast H_(∞) Filtering Algorithm

The calculation of the Riccati equation (covariance equation of a stateestimation error) shown in equation (14) is dominant in thecomputational complexity of the modified H_(∞) filter. Therefore, toprocess the modified H_(∞) filter at a high speed, if the filter gain ofequation (13) is directly determined without using the Riccati equation,the computational burden can be significantly reduced.

Since it is difficult to derive a fast algorithm for directly obtaininga filter gain K_(s,k)εR^(N×1), however, evolving an algorithm for fastcalculating a gain matrix defined as follows is examined.K_(k)=P_(k)C_(k) ^(T)εR^(N×2)  (16)where,

$\begin{matrix}{{p_{k} = {\left\lbrack {O_{k}^{T}\Omega_{k}O_{k}} \right\rbrack^{- 1} = \left\lbrack {\sum\limits_{i = 1}^{k}{\rho^{k - i}C_{i}^{T}W_{i}C_{i}}} \right\rbrack^{- 1}}}{{\Omega_{k} = \begin{bmatrix}{\rho\;\Omega_{k - 1}} & 0 \\0 & W_{k}\end{bmatrix}},{\Omega_{1} = W_{1}},{W_{i} = {{\rho\; R_{i}^{- 1}} = {\begin{bmatrix}1 & 0 \\0 & {- \gamma_{f}^{- 2}}\end{bmatrix} \in R^{2 \times 2}}}}}{{O_{k} = \begin{bmatrix}C_{1} \\\vdots \\C_{k}\end{bmatrix}},{C_{i} = {\begin{bmatrix}H_{i} \\H_{i}\end{bmatrix} \in {R^{2 \times N}.}}}}} & (17)\end{matrix}$

Here, the following lemmas are formed.

Lemma 1

A matrix P_(k) satisfies the Riccati equation of (14). Therefore, when again matrix K_(k) is obtained, the filter gain K_(s,k) is immediatelyobtained from the following lemma.

Lemma 2

The filter gain K_(s,k) of the modified H_(∞) filter is obtained byusing the gain matrix K_(k) as shown below. In practice, the gain matrixK_(k) can be fast calculated by the recursive method in Lemma 3.K _(s,k) =G _(k) ⁻¹ {tilde over (K)} _(k) , G _(k)=ρ+γ_(f) ⁻² H _(k){tilde over (K)} _(k) εR  (18)where,{tilde over (K)} _(k)(i)=ρK _(k)(i, 1), i=1, 2, . . . , N.  (19)Lemma 3

The gain matrix K_(k) is updated as follows.K _(k+1) =m _(k) −B _(k) F _(k) ⁻¹μ_(k) εR ^(N×2)  (20)

Here, m_(k)εR^(N×2) and μ_(k)εR^(1×2) are obtained by dividing a matrixof K^(U) _(k)=Q^(U) _(k) ⁻¹C^(U) _(k) as shown below.

$\begin{matrix}{\begin{bmatrix}m_{k} \\\mu_{k}\end{bmatrix} = {\begin{bmatrix}0 \\K_{k}\end{bmatrix} + {\begin{bmatrix}S_{k}^{- 1} \\{A_{k}S_{k}^{- 1}}\end{bmatrix}\left\lbrack {c_{k}^{T} + {A_{k}^{T}C_{k}^{T}}} \right\rbrack}}} & (21)\end{matrix}$

Auxiliary variables A_(k)εR^(N×1), S_(k)εR, and B_(k)F⁻¹ _(k)εR^(N×1)are obtained as well.

In conclusion, the fast H_(∞) filtering algorithm can be summarized asbelow.

FIG. 1 shows a flowchart of the fast algorithm, where L indicates amaximum data length.

[Step 0] Set initial conditions of a recursive expression as follows,where ε₀ is a substantially large positive constant.

${K_{0} = 0},{A_{- 1} = 0},{S_{- 1} = \frac{\rho}{ɛ_{0}}},{D_{- 1} = 0},{{\hat{x}}_{0|0} = 0}$[Step 1] Compare time k with the maximum data length L. When the time kis larger than the maximum data length, terminate the processing. Whenthe time k is equal to or smaller than the maximum data length, theprocessing proceeds to the next step (a conditional statement can beremoved, if unnecessary).[Step 2] Determine A_(k) and S_(k) recursively as follows.

$\begin{matrix}{{\overset{\sim}{e}}_{k} = {c_{k} + {C_{k}A_{k - 1}}}} & {\varepsilon\; R^{2 \times 1}} \\{A_{k} = {A_{k - 1} - {K_{k}W_{k}{\overset{\sim}{e}}_{k}}}} & {\varepsilon\; R^{N \times 1}} \\{e_{k} = {c_{k} + {C_{k}A_{k}}}} & {\varepsilon\; R^{2 \times 1}} \\{S_{k} = {{\rho\; S_{k - 1}} + {e_{k}^{T}W_{k}{\overset{\sim}{e}}_{k}}}} & {\varepsilon\; R}\end{matrix}\quad$[Step 3] Calculate K^(U) _(k) as follows.

${\overset{\Cup}{K}}_{k} = {\begin{bmatrix}{S_{k}^{- 1}e_{k}^{T}} \\{K_{k} + {A_{k}S_{k}^{- 1}e_{k}^{T}}}\end{bmatrix} \in R^{{({N + 1})} \times 2}}$[Step 4] Divide K^(U) _(k) as follows.

${{\overset{\Cup}{K}}_{k} = {{\begin{bmatrix}m_{k} \\\mu_{k}\end{bmatrix}\mspace{14mu} m_{k}} \in R^{N \times 2}}},{\mu_{k} \in R^{1 \times 2}}$[Step 5] Determine D_(k), and obtain a gain matrix K_(s,k+1) fromK_(k+1) as follows:η_(k) =c _(k−N) +C _(k+1) D _(k−1)D _(k) =[D _(k−1) −m _(k) W _(k)η_(k)][1−μ_(k) W _(k)η_(k)]⁻¹K _(k+1) =m _(k) −D _(k)μ_(k){tilde over (K)} _(k+1)(i)=ρK _(k+1)(i, 1), i=1, . . . , NK _(s,k+1) =G _(k+1) ⁻¹ {tilde over (K)} _(k+1) , G _(k+1)=ρ+γ_(f) ⁻² H_(k+1) {tilde over (K)} _(k+1)where, η_(k)εR^(2×1), D_(k)εR^(N×1), K_(k+1)εR^(N×2), K_(s,k+1)εR^(N×1),0<ρ=1−γ_(f) ⁻²≦1, γ_(f)>1.[Step 6] Update the filter equation of the H_(∞) filter as follows.{circumflex over (x)} _(k+1|k+1) ={circumflex over (x)} _(k|k) +K_(s,k+1)(y _(k+1) −H _(k+1) {circumflex over (x)} _(k|k))[Step 7] Put the time k forward (k=k+1). The processing returns to Step2, and the subsequent processes are repeated as long as the data exists.Lemma 6 (Existence Condition Suitable for Fast Processing)

The existence of the fast H_(∞) filter can be checked with thecomputational complexity of O(N) by using the following existencecondition.

[Existence Condition]−l{circumflex over (Ξ)}i+ργ _(f) ²>0, i=0, . . . , k  (22)where,

$\begin{matrix}{{\varrho = {1 - \gamma_{f}^{2}}},{{\hat{\Xi}}_{i} = \frac{H_{i}{\overset{\sim}{K}}_{i}}{1 - {H_{i}{\overset{\sim}{K}}_{i}}}}} & (23)\end{matrix}$4. Computational Complexity of the Present Fast Algorithm

Next, how the computational complexity of the fast H_(∞) filteringalgorithm decreases, as compared with the computational requirement ofthe modified H_(∞) filtering algorithm, will be discussed. Only thenumber of multiplications is used for evaluating the amount ofcalculation of an equation, and is calculated by the following method.Number of multiplications when a J-by-K matrix is multiplied by a K-by-Lmatrix is J×K×L (times).

Here, when three or more matrixes or vectors are multiplied, they arecalculated from the left unless a direction is specially shown in thefigure.

(Computational Complexity of the Modified H_(∞) Filtering Algorithm)

FIGS. 2 and 3 are views showing of the amount of calculation of eachpart of the modified H_(∞) filtering algorithm, where N indicates a tapnumber (the dimension of the state vector X_(k)). In FIG. 3( a), acalculation for obtaining R_(e,k) ⁻¹ from R_(e,k) is ignored. Similarly,in FIG. 2( a), a calculation for obtaining (H_(k+1)P^_(k+1|k)H^(T)_(k+1)+1)⁻¹ from (H_(k+1)P^_(k+1|k)H^(T) _(k+1)+1) is also ignored.

As shown in FIGS. 2( a), 3(a), and 3(b), the amount of calculation ofeach of K_(s|k+1), R_(e,k), and P^_(k+1|k) is in proportion to thesquare of the tap number. Therefore, the amount of calculation of theentire modified H_(∞) filtering algorithm is O(N²) per unit time step.

FIG. 4 is a view showing the amount of calculation required when theorder of matrix calculations is changed. More specifically, FIG. 4 showsthe amount of calculation required when the order of matrix calculationsof the second term on the right-hand side is changed in the Riccatiequation.

Since the amount of calculation of the above-described part isproportional to the cube of the tap number, the amount of calculation ofP^_(k+1|k) is also in proportion to the cube of the tap number.Accordingly, the amount of calculation of the entire H_(∞) filterincreases from the square to the cube of the tap number.

Since either algorithm requires the amount of calculation proportionalto the square or cube of the tap number, however, the computationalburden for carrying out the filter increases significantly as the tapnumber increases. In fact, since a tap number used in the field ofcommunication engineering, for example, is approximately 400, thepractical use of the algorithm becomes very difficult.

(Computational Complexity of the Fast H_(∞) Filtering Algorithm)

FIGS. 5 and 6 are views showing the amount of calculation in the fastH_(∞) filtering algorithm. In the expression of K^(u) _(k) in FIG. 5(b), S_(k) ⁻¹ is obtained from S_(k), but the calculation thereof isignored. Similarly, in the expression of D_(k) in FIG. 6( b), acalculation for obtaining [1−μ_(k)W_(k)η_(k)]⁻¹ from [1−μ_(k)W_(k)η_(k)]is also ignored.

The amount of calculation in the entire present fast algorithm is O(N)per unit time step according to FIGS. 5 and 6. Therefore, the amount ofcalculation in the fast H_(∞) filtering algorithm is in proportion tothe tap number. In this case, the amount of calculation (the number ofmultiplications) for performing the fast H_(∞) filter once is 28N+16 perunit step, and is approximately double the amount (multiplicationfrequency) of calculation required for a fast Kalman filter, that is12N+3.

As described above, although the computational complexity proportionalto the square or cube of the tap number is required for the modifiedH_(∞) filtering algorithm, the computational complexity of the presentfast algorithm is smaller and proportional to the tap number.

5. Echo Canceller

The advantage of the present invention will be examined, with an echocanceller being taken as an example.

An observation value {y_(k)} of an echo {d_(k)} is expressed in thefollowing expression by an (time-varying) impulse response {h_(i)[k]} ofan echo path, where it is considered that a received signal {u_(k)} isan input signal to the echo path:

$\begin{matrix}{{y_{k} = {{d_{k} + \upsilon_{k}} = {{\sum\limits_{i = 0}^{N - 1}{{h_{i}\lbrack k\rbrack}u_{k - i}}} + \upsilon_{k}}}},{k = 0},1,2,\ldots} & (24)\end{matrix}$where, u_(k) and y_(k) indicate, respectively, the received signal andthe echo at time t_(k) (=kT, T is a sampling period); v_(k) indicatescircuit noise having zero mean at time t_(k); and h_(i)[k] (i=0, . . . ,and N−1) is a time-varying impulse responses assuming a gradual change,and the tap number N thereof is known. Once estimated values {h^_(i)[k]}of the impulse response are obtained, a quasi echo is generated asfollows by using the estimated values.

$\begin{matrix}{{{\hat{d}}_{k} = {\sum\limits_{i = 0}^{N - 1}{{{\hat{h}}_{i}\lbrack k\rbrack}u_{k - i}}}},{k = 0},1,2,\ldots} & (25)\end{matrix}$

Subtracting this from the echo (y_(k)−d^_(k)≈0), the echo is cancelled,where u_(k−1)=0 when k−i<0.

From the above description, the echo canceller problem is equivalent tosuccessively estimating the impulse response {h_(i)[k]} of the echo pathfrom the received signal {u_(k)} and echo {y_(k)}, both of which aredirectly observable.

In general, when the H_(∞) filter is applied to an echo canceller,equation (24) has to be expressed by a state-space model formed of astate equation and an observation equation. In this case, since thestate vector to be obtained is the impulse response {h_(i)[k]}, allowinga state vector x_(k) to fluctuate with w_(k), the following state-spacemodel can be constructed for the echo path.

$\begin{matrix}{{x_{k + 1} = {x_{k} + w_{k}}},} \\{{y_{k} = {{H_{k}x_{k}} + v_{k}}},} \\{{z_{k} = {H_{k}x_{k}}},}\end{matrix}\mspace{140mu}\begin{matrix}{x_{k},{w_{k}\varepsilon\; R^{N}}} \\{y_{k},{v_{k}\varepsilon\; R}} \\{{z_{k}\varepsilon\; R},{H_{k}\varepsilon\; R^{1 \times N}}}\end{matrix}\mspace{155mu}\begin{matrix}{\text{(26)}} \\{\text{(27)}} \\{\text{(28)}}\end{matrix}$where,x _(k) =[h ₀ [k], . . . , h _(N−1) [k]] ^(T) , w _(k) =[w _(k)(1), . . ., w _(k)(N)]^(T)H _(k) =[u _(k) , . . . , u _(k−N+1)], L_(k)=H_(k)

Modified H_(∞) filtering algorithm and fast H_(∞) filtering algorithmfor such a state-space model are the same as those described above.While the impulse response is estimated, if the occurrence of atransmission signal is detected, the estimation is generally stopped inthe meanwhile.

Thus, when an estimate {h^_(i)[k]} of the impulse response is obtained,the quasi echo is successively obtained therefrom as follows.

$\begin{matrix}{{\hat{d}}_{k} = {{H_{k}{\hat{x}}_{k❘k}} = {\sum\limits_{i = 0}^{N - 1}{{{\hat{h}}_{i}\lbrack k\rbrack}u_{k - i}}}}} & (29)\end{matrix}$Therefore, subtracting this from an actual echo to cancel the echo, anecho canceller is implemented. Here, an estimate error,e_(k)=y_(k)−d^_(k), is called a residual echo.6. Evaluation for Time-Invariant Impulse Response(Evaluation of Estimation Accuracy)

A modified H_(∞) filter and a fast H_(∞) filter are evaluated bysimulation in a case in which the impulse response of an echo path istime-invariant (h_(i)[k]=h_(i)) and the tap number N thereof is 24.

$\begin{matrix}{y_{k} = {{\sum\limits_{i = 0}^{23}{h_{i}u_{k - i}}} + \upsilon_{k}}} & (30)\end{matrix}$

FIG. 7 is a view showing values of the impulse response {h_(i)} in thiscase. v_(k) is stationary Gaussian white noise having zero mean andvariance σ_(v) ² of 1.0×10⁻⁶, and a sampling period T is set to 1.0 forconvenience.

The received signal {u_(k)} is approximated by a quadratic AR model asshown below.u _(k)=α₁ u _(k−1)+α₂ u _(k−2) +w _(k)′  (31)where, α₁=0.7, α₂=0.1, and w_(k)′ is stationary Gaussian white noisehaving zero mean and variance σ_(w′) ² of 0.04.

The modified H_(∞) filter and the fast H_(∞) filter will be compared.

FIG. 8 includes views showing estimated results of the impulse responsesof the modified H_(∞) filter and the fast H_(∞) filter (initial valuex^_(0|0)=0, ε_(o)=20). FIGS. 8( a) and (b) show estimated results ofboth filters when γ_(f)=10⁵, and FIGS. 8( c) and (d) show estimatedresults (x^_(100|100)) thereof when γ_(f)=2.0. From the figures,performance on the estimation accuracy of both filters is equal. Inother words, speeding-up does not reduce the estimation accuracy. Notethat, if γ_(f) is too small, the existence condition of the filters isnot satisfied. When γ_(f)=1.0×10⁵, the results are substantially equalto that of a fast Kalman filter. Therefore, it is found that the fastH_(∞) filtering algorithm includes the fast Kalman filtering algorithmand its convergence rate can be accelerated by adjusting γ_(f).

(Evaluation of Computation Time)

Next, the computation time required for the modified H_(∞) filter andthat for the fast H_(∞) filter are evaluated under conditions where theimpulse response of the echo path is time-invariant and the tap numberis increased to 24, 48, 96, 192, and 384. Since dispersion may occur inone measurement, the average of four measurements was used. The valuesshown in FIG. 7 are used as impulse responses {h_(i)} in simulation, andimpulse responses {h_(i)} thereafter (24≦k<N) are set to 0. The filtercalculation is performed up to step 100. The computation time wasmeasured by a command “etime” of MATLAB on a workstation (sparc, 60 MHz,32 MB).

FIG. 9 is a view showing measurement results of the computation time. Inthe Riccati equation, matrix calculation is performed for a modifiedH_(∞) filter (2) such that the amount of calculation is in proportion tothe square of the tap number, and matrix calculation is performed for amodified H_(∞) filter (1) such that the amount of calculation is inproportion to the cube of the tap number (see FIG. 3( b) and FIG. 4). Inmodified H_(∞) filters, since the computational complexity is inproportion to the square or cube of the tap number depending on theorder of matrix calculation as described above, they are not practical.

7. Evaluation for Time-Varying Impulse Response

(Evaluation of Tracking Performance)

The tracking performance of each algorithm will be evaluated by usingthe echo canceller in a case in which the system (impulse response) isvaried with time. It is assumed that the tap number of the impulseresponse is 48, and {h_(i)} is varied with time, as shown in FIG. 10(a), based on the values shown in FIG. 7. It is also assumed that v_(k)is stationary Gaussian white noise having zero mean and variance σ_(v) ²of 1.0×10⁻⁶, and the sampling period T is set to 1 for convenience. Thereceived signal {u_(k)} is approximated by a quadratic AR model asfollows.u _(k)=α₁ u _(k−1)+α₂ u _(k−2) +w _(k)′  (32)

Here, α₁=0.7, α₂=0.1, and w_(k)′ is stationary Gaussian white noisehaving zero mean and variance σ² _(w′) of 0.04.

FIGS. 10 and 11 are views showing the simulation result of eachalgorithm. These views show the tracking performance of time-varyingsystems which employ a fast H_(∞) filter (fast HF), a fast Kalman filter(fast KF), and LMS algorithm (LMS). FIG. 10( b) shows the estimatesobtained with the fast H_(∞) filter when γ_(f)=2.0. FIG. 11( a) showsthe estimates obtained with the fast Kalman filter. The initial value ofthe fast H_(∞) filter is set such that x^_(0|0)=0 and ε_(o)=20, and theinitial value of the fast Kalman filter is set in the same way. FIG. 11(b) shows the estimates obtained by the LMS algorithm, wherein theinitial value is set such that h^_(o)=0, and the step size is set suchthat μ=0.5 so as to give a stable and rapid convergence. It is foundthat the tracking performance of the fast H_(∞) filter is extremelyexcellent, and the estimates become stable in about thirty steps afterthe impulse response is varied. On the other hand, the fast Kalmanfilter and the LMS algorithm cannot track the impulse response at all.

Generally, the tracking performance of H_(∞) filters having no systemnoise drops with time since the filter gain becomes smaller due to adecay in the diagonal component of P^_(k|k−1) and the amount of updateof the estimates decreases. In other words, as the number of stepsincreases, the estimates are updated little. Therefore, in order toimprove the tracking performance of Kalman filters and H_(∞) filters, anappropriate value needs to be externally added to the diagonal componentof the matrix P^_(k|k−1). If it is directly added, however, a fastalgorithm which uses the shift property of an observation matrix H_(k)cannot be implemented. It is one of significant features of the presentinvention to solve this problem theoretically by applying a weight ρ of1−γ_(f) ⁻² to the H_(∞) evaluation criterion. The weight ρ appears in anupdate equation of S_(k) of the fast H_(∞) filtering algorithm, asfollows.

(Update of Auxiliary Variable S_(k) of the Fast H_(∞) Filter)

An auxiliary variable S_(k) of the fast H_(∞) filter is indicated by thefollowing expression.S _(k) =ρS _(k−1) +e ^(T) _(k) W _(k) e{tilde over ( )} _(k),0<ρ=1−γ_(f) ⁻²≦1

In the fast H_(∞) filtering algorithm, S_(k) is used as S_(k) ⁻¹ in theequation of K^(u) _(k). In order to largely update the filter equation,S_(k) ⁻¹ must be larger. In other word, S_(k) needs to be kept small tomake the large update. The existence of ρ prevents S_(k) from increasingrapidly, which is resultantly equivalent to adding system noise, andthereby the tracking performance is improved. Since the weight ρ isdefined as 1−γ_(f) ⁻², the tracking performance can be varied byadjusting γ_(f) as confirmed in the simulation.

FIG. 12 is a view showing the relationship between γ_(f) and ρ.According to the figure, when γ_(f)=3.0, ρ=0.8889, which means that 89%of S_(k−1) is transmitted to S_(k). Note that, if γ_(f) is set verysmall, however, the effect of S_(k−1) is significantly reduced and theexistence condition of the filter is not satisfied. When γ_(f) is large,ρ1. An increase in S_(k) is not suppressed at all, and therefore, thetracking performance drops. When γ_(f)=∞(ρ=1, in particular, the presentfast algorithm completely matches the fast Kalman filtering algorithm.

(Evaluation of Computation Time)

FIG. 13 is a view showing the relationship among the tap number and thecomputation time for the fast H_(∞) filter, the fast Kalman filter, andthe LMS algorithm, where the number of time steps executed for thefilters is 300 and γ_(f)=3.0. The computation time was measured for thefast H_(∞) filter, the fast Kalman filtering algorithm, and the LMSalgorithm when the tap number was increased to 48, 96, 192, and 384 inthe cases shown in FIGS. 10 and 11. Because dispersion may occur in onemeasurement result, the average of four measurement results, [forexample,] was used.

It can be confirmed that, in any algorithm, the amount of calculation isin proportion to the tap number. It is also found that when the tapnumber is large, the computation time for the fast H_(∞) filteringalgorithm is about a little less than twice the computation time for thefast Kalman filtering algorithm, and is approximately four times longerthan that for the LMS algorithm, which is practical. Considering thetracking performance, it can be said that the fast H_(∞) filteringalgorithm is sufficiently effective.

8. Demonstration of Lemmas

Now, the above-described lemmas will be demonstrated.

(Demonstration of Lemma 1)

The inverse matrix of P_(k) will be indicated by equation (33). Further,a recursive equation for the matrix P_(k) can be obtained, as shown inequation (34), by using the matrix inversion lemma.

$\begin{matrix}\begin{matrix}{P_{k}^{- 1} = {{\rho\; O_{k - 1}^{T}\Omega_{k - 1}O_{k - 1}} + {C_{k}^{T}W_{k}C_{k}}}} \\{= {{\rho\; P_{k - 1}^{- 1}} + {C_{k}^{T}W_{k}{C_{k}.}}}}\end{matrix} & (33) \\\begin{matrix}{P_{k} = \left\lbrack {{\rho\; P_{k - 1}^{- 1}} + {\left\lbrack {H_{k}^{T}\mspace{14mu} H_{k}^{T}} \right\rbrack{W_{k}\begin{bmatrix}H_{k} \\H_{k}\end{bmatrix}}}} \right\rbrack^{- 1}} \\{= {{\rho^{- 1}P_{k - 1}} - {\rho^{- 1}{{P_{k - 1}\left\lbrack {H_{k}^{T}\mspace{14mu} H_{k}^{T}} \right\rbrack} \cdot}}}} \\{{{\left( {W_{k}^{- 1} + {\begin{bmatrix}H_{k} \\H_{k}\end{bmatrix}\rho^{- 1}{P_{k - 1}\left\lbrack {H_{k}^{T}\mspace{14mu} H_{k}^{T}} \right\rbrack}}} \right)^{- 1} \cdot \begin{bmatrix}H_{k} \\H_{k}\end{bmatrix}}\rho^{- 1}P_{k - 1}},} \\{{\rho\; P_{k}} = {P_{k - 1} - {{P_{k - 1}\left\lbrack {H_{k}^{T}\mspace{14mu} H_{k}^{T}} \right\rbrack} \cdot}}} \\{{{\left( {R_{k} + {\begin{bmatrix}H_{k} \\H_{k}\end{bmatrix}{P_{k - 1}\left\lbrack {H_{k}^{T}\mspace{14mu} H_{k}^{T}} \right\rbrack}}} \right)^{- 1} \cdot \begin{bmatrix}H_{k} \\H_{k}\end{bmatrix}}P_{k - 1}},} \\{P_{k} = {P_{k - 1} - {{P_{k - 1}\left\lbrack {H_{k}^{T}\mspace{14mu} H_{k}^{T}} \right\rbrack} \cdot}}} \\{{{\left( {R_{k} + {\begin{bmatrix}H_{k} \\H_{k}\end{bmatrix}{P_{k - 1}\left\lbrack {H_{k}^{T}\mspace{14mu} H_{k}^{T}} \right\rbrack}}} \right)^{- 1} \cdot \begin{bmatrix}H_{k} \\H_{k}\end{bmatrix}}P_{k - 1}} + {\gamma_{f}^{- 2}{P_{k}.}}}\end{matrix} & (34)\end{matrix}$

It is understood, when P_(k) is replaced with P^_(k+1|k), that the aboveequation satisfies the Riccati equation of (14).

(Demonstration of Lemma 2)

The gain matrix K_(k) can be expressed as follows.

$\begin{matrix}\begin{matrix}{K_{k} = {{P_{k}C_{k}^{T}} = {\left\lbrack {{\rho\; P_{k - 1}^{- 1}} + {C_{k}^{T}W_{k}C_{k}}} \right\rbrack^{- 1}C_{k}^{T}}}} \\{= {{\rho^{- 1}P_{k - 1}C_{k}^{T}} - {\rho^{- 1}P_{k - 1}{C_{k}^{T} \cdot \left\lbrack {W_{k}^{- 1} + {C_{k}\rho^{- 1}P_{k - 1}C_{k}^{T}}} \right\rbrack^{- 1}}}}} \\{C_{k}\rho^{- 1}P_{k - 1}C_{k}^{T}} \\{= {{\rho^{- 1}P_{k - 1}C_{k}^{T}} - {\rho^{- 1}P_{k - 1}{{C_{k}^{T}\left\lbrack {W_{k}^{- 1} + {C_{k}\rho^{- 1}P_{k - 1}C_{k}^{T}}} \right\rbrack}^{- 1} \cdot}}}} \\{\left\lbrack {\left( {W_{k}^{- 1} + {C_{k}\rho^{- 1}P_{k - 1}C_{k}^{T}}} \right) - W_{k}^{- 1}} \right\rbrack} \\{= {\rho^{- 1}P_{k - 1}{C_{k}^{T}\left\lbrack {I + {W_{k}C_{k}\rho^{- 1}P_{k - 1}C_{k}^{T}}} \right\rbrack}^{- 1}}} \\{= {\rho^{- 1}P_{k - 1}C_{k}^{T}{W_{k} \cdot \left\lbrack {W_{k} + {\rho^{- 1}W_{k}C_{k}P_{k - 1}C_{k}^{T}W_{k}}} \right\rbrack^{- 1}}}} \\{= {\rho^{- 1}{{P_{k - 1}\left\lbrack {H_{k}^{T}\mspace{14mu} - {\gamma_{f}^{- 2}H_{k}^{T}}} \right\rbrack}\left\lbrack {\begin{bmatrix}1 & 0 \\0 & {- \gamma_{f}^{- 2}}\end{bmatrix} +} \right.}}} \\\left. {{\rho^{- 1}\begin{bmatrix}H_{k} \\{{- \gamma_{f}^{- 2}}H_{k}}\end{bmatrix}}{P_{k - 1}\left\lbrack {H_{k}^{T}\mspace{14mu} - {\gamma_{f}^{- 2}H_{k}^{T}}} \right\rbrack}} \right\rbrack^{- 1} \\{= {\rho^{- 1}{P_{k - 1}\left\lbrack {H_{k}^{T}\mspace{14mu} H_{k}^{T}} \right\rbrack}\left( {1 + {H_{k}P_{k - 1}H_{k}^{T}}} \right)^{- 1}}}\end{matrix} & (35)\end{matrix}$

Further, the filter gain can be obtained from the first block column ofthe gain matrix K_(k), as shown in equation (18), by usingG_(k)=(ρ+H_(k)P_(k−1)H^(T) _(k))/(1+H_(k)P_(k−1)H^(T) _(k)) andH_(k)K{tilde over ( )}_(k)=H_(k)P_(k−1)H^(T) _(k)/(1+H_(k)P_(k−1)H^(T)_(k)).

(Demonstration of Lemma 3)

Assuming that the gain matrix K_(i) (i=0, . . . , and k) is given, thefollowing matrix, K_(k+1), will be calculated.Q _(k+1) K _(k+1) =C _(k+1) ^(T)  (36)First, equations (37) and (38) are newly introduced to utilize the shiftproperty of C_(k). Q^(U) _(k) is expressed recursively as shown inequation (39), and is divided as in the following equation (40).

$\begin{matrix}{{\overset{\Cup}{C}}_{k}^{T} = {\begin{bmatrix}c_{k}^{T} \\C_{k}^{T}\end{bmatrix} = {\begin{bmatrix}C_{k + 1}^{T} \\c_{k - N}^{T}\end{bmatrix} \in R^{{({N + 1})} \times 2}}}} & (37) \\{{\overset{\Cup}{Q}}_{k} = {{\sum\limits_{i = 1}^{k}{\rho^{k - i}{\overset{\Cup}{C}}_{i}^{T}W_{i}{\overset{\Cup}{C}}_{i}}} \in R^{{({N + 1})} \times {({N + 1})}}}} & (38) \\{{\overset{\Cup}{Q}}_{k} = {{\rho\;{\overset{\Cup}{Q}}_{k - 1}} + {{\overset{\Cup}{C}}_{k}^{T}W_{k}{{\overset{\Cup}{C}}_{k}.}}}} & (39) \\{{\overset{\Cup}{Q}}_{k} = {\begin{bmatrix}M_{k} & T_{k}^{T} \\T_{k} & Q_{k}\end{bmatrix} = {\begin{bmatrix}Q_{k + 1} & {\underset{—}{T}}_{k}^{T} \\{\underset{—}{T}}_{k} & {\underset{—}{M}}_{k}\end{bmatrix}.}}} & (40)\end{matrix}$

Using this notation, equation (36) of the time steps k and k+1 isincluded in the following equation.

$\begin{matrix}{{{\overset{\Cup}{Q}}_{k}\begin{bmatrix}0 \\K_{k}\end{bmatrix}} = {\begin{bmatrix}\alpha_{k}^{T} \\C_{k}^{T}\end{bmatrix} = {{\overset{\Cup}{C}}_{k}^{T} + \begin{bmatrix}{\alpha_{k}^{T} - c_{k}^{T}} \\0\end{bmatrix}}}} & (41) \\{{{\overset{\Cup}{Q}}_{k}\begin{bmatrix}K_{k + 1} \\0\end{bmatrix}} = {\begin{bmatrix}C_{k + 1}^{T} \\\beta_{k}^{T}\end{bmatrix} = {{\overset{\Cup}{C}}_{k}^{T} + \begin{bmatrix}0 \\{\beta_{k}^{T} - c_{k - N}^{T}}\end{bmatrix}}}} & (42)\end{matrix}$

-   -   where,        α_(k) ^(T) =T _(k) ^(T) K _(k) εR ^(1×2), β_(k) ^(T) =T _(k) K        _(k+1) εR ^(1×2.)

Based on the notation, it is more convenient to obtain K^(U)_(k)εR^((N+1)×2), which satisfies the following equation, than to obtainK_(k) directly.{hacek over (Q)}_(k){hacek over (K)}_(k)={hacek over (C)}_(k) ^(T)  (43)where,{hacek over (K)} _(k) =[k _(k+1) ^(T) K _(k) ^(T)]^(T) =[K _(k+1) ^(T) k_(k−N) ^(T)]^(T).  (44)

To this end, K^(U) _(k)εR^((N+1)×2) can be expressed as shown inequation (46) by using equation (45), obtained from equation (41).

$\begin{matrix}{{\overset{\Cup}{C}}_{k}^{T} = {{{\overset{\Cup}{Q}}_{k}\begin{bmatrix}0 \\K_{k}\end{bmatrix}} - \begin{bmatrix}{\alpha_{k}^{T} - c_{k}^{T}} \\0\end{bmatrix}}} & (45) \\\begin{matrix}{{\overset{\Cup}{K}}_{k} = {\begin{bmatrix}m_{k} \\\mu_{k}\end{bmatrix} = {{{\overset{\Cup}{Q}}_{k}^{- 1}{\overset{\Cup}{C}}_{k}^{T}} = {\begin{bmatrix}0 \\K_{k}\end{bmatrix} - {{\overset{\Cup}{Q}}_{k}^{- 1}\begin{bmatrix}{\alpha_{k}^{T} - c_{k}^{T}} \\0\end{bmatrix}}}}}} \\{= {\begin{bmatrix}0 \\K_{k}\end{bmatrix} - {\begin{bmatrix}S_{k}^{- 1} \\{A_{k}S_{k}^{- 1}}\end{bmatrix}\left\lbrack {\alpha_{k}^{T} - c_{k}^{T}} \right\rbrack}}}\end{matrix} & (46)\end{matrix}$

Here, K^(U) _(k) is divided into m_(k)εR^(N×2) and μ_(k)εR^(1×2). Alsonote that α^(T) _(k)−c^(T) _(k)=−(c_(k) ^(T)+A_(k) ^(T)C_(k) ^(T)).Further, assuming that Q^(U) _(k) has an inverse matrix, auxiliaryvariables A_(k)εR^(N×1) and S_(k)εR satisfy the following equation.

$\begin{matrix}{{{\overset{\Cup}{Q}}_{k}\begin{bmatrix}1 \\A_{k}\end{bmatrix}} = {\begin{bmatrix}S_{k} \\0\end{bmatrix}\left( {{\begin{bmatrix}1 \\A_{k}\end{bmatrix}S_{k}^{- 1}} = {{\overset{\Cup}{Q}}_{k}^{- 1}\begin{bmatrix}1 \\0\end{bmatrix}}} \right)}} & (47)\end{matrix}$where, the bottom block of the above equation means T_(k)+Q_(k)A_(k)=0or T_(k) ^(T)=−A_(k) ^(T)Q_(k) ^(T).

Next, auxiliary variables B_(k)εR^(N×1) and F_(k)εR such as those shownin the following equation (48) are introduced to delete μ_(k) inequation (46) without affecting the top block of C^(T) _(k). Further,subtracting B^(U) _(k)F_(k) ⁻¹μ_(k) from K^(U) _(k) in equation (46)provides equation (49).

$\begin{matrix}{{{\overset{\Cup}{Q}}_{k}{\overset{\Cup}{B}}_{k}} = {{{\overset{\Cup}{Q}}_{k}\begin{bmatrix}B_{k} \\F_{k}\end{bmatrix}} = {\begin{bmatrix}0 \\1\end{bmatrix}\left( {{\overset{\Cup}{B}}_{k} = \begin{bmatrix}B_{k} \\F_{k}\end{bmatrix}} \right)}}} & (48) \\\begin{matrix}{{{\overset{\Cup}{K}}_{k} - {{\overset{\Cup}{B}}_{k}F_{k}^{- 1}\mu_{k}}} = {\begin{bmatrix}m_{k} \\\mu_{k}\end{bmatrix} - {\begin{bmatrix}{B_{k}F_{k}^{- 1}} \\1\end{bmatrix}\mu_{k}}}} \\{= \begin{bmatrix}{m_{k} - {B_{k}F_{k}^{- 1}\mu_{k}}} \\0\end{bmatrix}}\end{matrix} & (49)\end{matrix}$

Then, the left-hand side of equation (49) is multiplied by Q^(U) _(k)from the left to obtain the following equation.

$\begin{matrix}\begin{matrix}{{{\overset{\Cup}{Q}}_{k}\left( {{\overset{\Cup}{K}}_{k} - {{\overset{\Cup}{B}}_{k}F_{k}^{- 1}\mu_{k}}} \right)} = {{{\overset{\Cup}{Q}}_{k}{\overset{\Cup}{K}}_{k}} - {{\overset{\Cup}{Q}}_{k}{\overset{\Cup}{B}}_{k}F_{k}^{- 1}\mu_{k}}}} \\{= {{\overset{\Cup}{C}}_{k}^{T} - {\begin{bmatrix}0 \\1\end{bmatrix}F_{k}^{- 1}\mu_{k}}}} \\{= {{\overset{\Cup}{C}}_{k}^{T} - \begin{bmatrix}0 \\{F_{k}^{- 1}\mu_{k}}\end{bmatrix}}}\end{matrix} & (50)\end{matrix}$

Equation (49) is substituted for the left-hand side of the aboveequation. Then, equation (43) is expressed as follows:

$\begin{matrix}{{{{\overset{\Cup}{Q}}_{k}\left( {{\overset{\Cup}{K}}_{k} - {{\overset{\Cup}{B}}_{k}F_{k}^{- 1}\mu_{k}}} \right)} = {{\overset{\Cup}{C}}_{k}^{T} - \begin{bmatrix}0 \\{F_{k}^{- 1}\mu_{k}}\end{bmatrix}}},{{\begin{bmatrix}Q_{k + 1} & {\underset{—}{T}}_{k}^{T} \\{\underset{—}{T}}_{k} & {\underset{—}{M}}_{k}\end{bmatrix}\begin{bmatrix}{m_{k} - {B_{k}F_{k}^{- 1}\mu_{k}}} \\0\end{bmatrix}} = {\begin{bmatrix}C_{k + 1}^{T} \\c_{k - N}^{T}\end{bmatrix} + \begin{bmatrix}0 \\{{- F_{k}^{- 1}}\mu_{k}}\end{bmatrix}}}} & (51)\end{matrix}$This is the same form as equation (42). The following equation (52) canbe obtained from the top block of equation (51).Q _(k+1)(m _(k) −B _(k) F _(k) ⁻¹μ_(k))=C _(k+1) ^(T)  (52)Equations (36) and (52) are compared to obtain the update equation ofthe gain matrix K_(k).(Lemma 4)

The auxiliary variables A_(k) and S_(k) can be obtained as follows:A _(k) =A _(k−1) −K _(k) W _(k) [c _(k) +C _(k) A _(k−1) ]εR^(N×1)  (53)S _(k) =ρS _(k−1) +[c _(k) ^(T) +A _(k) ^(T) C _(k) ^(T) ]W _(k) [c _(k)+C _(k) A _(k−1) ]εR  (54)where, A⁻¹=0, and S⁻¹=1/ε₀.(Demonstration) By using the equation (55) of A_(k) and S_(k) andequation (39), equation (56) is obtained.

$\begin{matrix}{{{\overset{\Cup}{Q}}_{k - 1}\begin{bmatrix}1 \\A_{k - 1}\end{bmatrix}} = \begin{bmatrix}S_{k - 1} \\0\end{bmatrix}} & (55) \\\begin{matrix}{{{\overset{\Cup}{Q}}_{k}\begin{bmatrix}1 \\A_{k - 1}\end{bmatrix}} = {{\rho{{\overset{\Cup}{Q}}_{k - 1}\begin{bmatrix}1 \\A_{k - 1}\end{bmatrix}}} + {{\overset{\Cup}{C}}_{k}^{T}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}}}} \\{= {\begin{bmatrix}{\rho\; S_{k - 1}} \\0\end{bmatrix} + {\begin{bmatrix}c_{k}^{T} \\C_{k}^{T}\end{bmatrix}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}}}}\end{matrix} & (56)\end{matrix}$

On the other hand, the following equation is obtained by multiplyingboth sides of equation (41) by W_(k).

$\begin{matrix}{{{{\overset{\Cup}{Q}}_{k}\begin{bmatrix}0 \\K_{k}\end{bmatrix}}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}} = {\begin{bmatrix}\alpha_{k} \\C_{k}^{T}\end{bmatrix}{{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}.}}} & (57)\end{matrix}$By subtracting equation (57) from equation (56), the following equation(58) is formed.

$\begin{matrix}{{{{\overset{˘}{Q}}_{k}\left\lbrack {\begin{bmatrix}1 \\A_{k - 1}\end{bmatrix} - {\begin{bmatrix}0 \\K_{k}\end{bmatrix}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}}} \right\rbrack} = \mspace{59mu}{\begin{bmatrix}{\rho\; S_{k - 1}} \\0\end{bmatrix} + {\begin{bmatrix}c_{k}^{T} \\C_{k}^{T}\end{bmatrix}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}} - {\begin{bmatrix}\alpha_{k}^{T} \\C_{k}^{T}\end{bmatrix}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}}}},{{{\overset{˘}{Q}}_{k}\left\lbrack \;\begin{matrix}1 \\{A_{k - 1} - {K_{k}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}}}\end{matrix} \right\rbrack} = \mspace{284mu}\begin{bmatrix}{{\rho\; S_{k - 1}} + {\left\lbrack {c_{k}^{T} - \alpha_{k}^{T}} \right\rbrack{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}}} \\0\end{bmatrix}}} & (58)\end{matrix}$This equation is compared with equation (47). Since α_(k) ^(T)=T_(k)^(T)K_(k)=−A_(k) ^(T)C_(k) ^(T), equations (53) and (54) are obtained.(Lemma 5)

The auxiliary variable D_(k)=B_(k)F_(k) ⁻¹ is obtained by the followingequation (59). F_(k) is updated by the following equation (60).D _(k) =[D _(k−1) −m _(k) W _(k)η_(k)][1−μ_(k) W _(k)η_(k)]⁻¹ εR^(N×1)  (59)F _(k) =F _(k−1)[1−μ_(k) W _(k)η_(k) ]/ρεR  (60)where, η_(k)=C^(U) _(k)D^(U) _(k−1)=C_(k−N)+C_(k+1)D_(k−1), D⁻¹=0, andF⁻¹=0.(Demonstration) In order to update B_(k) and F_(k), equation (62) isformed by using equation (61).

$\begin{matrix}{{{\overset{˘}{Q}}_{k - 1}{\overset{˘}{B}}_{k - 1}} = {{{\overset{˘}{Q}}_{k - 1}\begin{bmatrix}B_{k - 1} \\F_{k - 1}\end{bmatrix}} = \begin{bmatrix}0 \\1\end{bmatrix}}} & (61) \\\begin{matrix}{{{\overset{˘}{Q}}_{k}{\overset{˘}{B}}_{k - 1}} = {{\rho{\overset{˘}{Q}}_{k - 1}{\overset{˘}{B}}_{k - 1}} + {{\overset{˘}{C}}_{k}^{T}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{B}}_{k - 1}}}} \\{= {{\rho\begin{bmatrix}0 \\1\end{bmatrix}} + {{\overset{˘}{C}}_{k}^{T}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{B}}_{k - 1}}}}\end{matrix} & (62)\end{matrix}$In order to modify the above equation so as to have the same form asequation (61), C^(U) _(k) ^(T)W_(k)C^(U) _(k)B_(k−1) ^(U) is subtractedfrom equation (62) to obtain the following equation.

$\begin{matrix}{{{{{\overset{˘}{Q}}_{k}{\overset{˘}{B}}_{k - 1}} - {{\overset{˘}{C}}_{k}^{T}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{B}}_{k - 1}}} = {{{{\overset{˘}{Q}}_{k}{\overset{˘}{B}}_{k - 1}} - {{\overset{˘}{Q}}_{k}{\overset{˘}{K}}_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{B}}_{k - 1}}} = {\rho\begin{bmatrix}0 \\1\end{bmatrix}}}},{{{\overset{˘}{Q}}_{k}\left\lbrack {{\overset{˘}{B}}_{k - 1} - {{\overset{˘}{K}}_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{B}}_{k - 1}}} \right\rbrack} = {\rho\begin{bmatrix}0 \\1\end{bmatrix}}}} & (63)\end{matrix}$Comparing the above last equation with equation (48) yields a recursiveequation for B^(U) _(k).{hacek over (B)} _(k)=({hacek over (B)} _(k−1) −{hacek over (K)} _(k) W_(k) {hacek over (C)} _(k) {hacek over (B)} _(k−1))/ρ  (64)

$\begin{matrix}{{D_{k} = {B_{k}F_{k}^{- 1}}},{{\overset{˘}{D}}_{k} = {{{\overset{˘}{B}}_{k}F_{k}^{- 1}} = \begin{bmatrix}D_{k} \\1\end{bmatrix}}}} & (65)\end{matrix}$B_(k) and F_(k) are updated by this equation.

Since they appear only as a combination of D^(U) _(k)=B^(U) _(k)F_(k) ⁻¹or D_(k)=B_(k)F_(k) ⁻¹ εR^(N×1), however, it is more convenient toexpress equations (48) and (64) with equation (65). The matrix D_(k)satisfies the following equation (66).

$\begin{matrix}{{{{\overset{˘}{Q}}_{k}{\overset{˘}{D}}_{k}} = {{{\overset{˘}{Q}}_{k}{\overset{˘}{B}}_{k}F_{k}^{- 1}} = {\begin{bmatrix}0 \\1\end{bmatrix}F_{k}^{- 1}}}},{{{\overset{˘}{Q}}_{k}\begin{bmatrix}D_{k} \\1\end{bmatrix}} = \begin{bmatrix}0 \\F_{k}^{- 1}\end{bmatrix}}} & (66) \\\begin{matrix}{\left. {{{\overset{˘}{Q}}_{k}{\overset{˘}{\left\lbrack B \right.}}_{k - 1}F_{k - 1}^{- 1}} - {{\overset{˘}{K}}_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{B}}_{k - 1}F_{k - 1}^{- 1}}} \right\rbrack = {{\overset{˘}{Q}}_{k}\left\lbrack {\overset{˘}{D}\;}_{k - 1} \right.}} \\\left. {{\overset{˘}{K}}_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{D}}_{k - 1}} \right\rbrack \\{= \begin{bmatrix}0 \\{\rho\; F_{k - 1}^{- 1}}\end{bmatrix}}\end{matrix} & (67)\end{matrix}$Next, equation (63) is multiplied by F_(k−1) ⁻¹ to obtain equation (67),and is further expressed by the following equation (68) when D^(U)_(k−1)=B^(U) _(k−1)F_(k−1) ⁻¹ is used.

$\begin{matrix}{{{{\overset{˘}{Q}}_{k}\left\lbrack {{\overset{˘}{D}}_{k - 1} - {\begin{bmatrix}m_{k} \\\mu_{k}\end{bmatrix}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{D}}_{k - 1}}} \right\rbrack} = \begin{bmatrix}0 \\{\rho\; F_{k - 1}^{- 1}}\end{bmatrix}},{{{\overset{˘}{Q}}_{k}\begin{bmatrix}{D_{k - 1} - {m_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{D}}_{k - 1}}} \\{1 - {\mu_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{D}}_{k - 1}}}\end{bmatrix}} = \begin{bmatrix}0 \\{\rho\; F_{k - 1}^{- 1}}\end{bmatrix}}} & (68)\end{matrix}$Therefore, the following equation is obtained when equation (68) ismultiplied by [1−μ_(k)W_(k)C^(U) _(k)D^(U) _(k−1)]⁻¹.

${{\overset{˘}{Q}}_{k}\begin{bmatrix}{\left\lbrack {D_{k - 1} - {m_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{D}}_{k - 1}}} \right\rbrack\left\lbrack {1 - {\mu_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{D}}_{k - 1}}} \right\rbrack}^{- 1} \\1\end{bmatrix}} = \mspace{464mu}\begin{bmatrix}0 \\{\rho\;{F_{k - 1}^{- 1}\left\lbrack {1 - {\mu_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{D}}_{k - 1}}} \right\rbrack}^{- 1}}\end{bmatrix}$By comparing this equation with equation (66), an update equation forD_(k) and F_(k) is finally obtained.(Demonstration of Lemma 6 (Existence Condition Suitable for FastProcessing))

As described above, the existence of the fast H_(∞) filter can bechecked with the computational complexity of O(N) by using the existencecondition of equations (22) and (23). A demonstration thereof will beshown below.

When the characteristic equation of a 2×2 matrix R_(e,k) shown in thefollowing equation (69) is solved, a eigenvalue λ_(i) of R_(e,k) isobtained by the following equation (70).

$\begin{matrix}\begin{matrix}{{{{\lambda\; I} - R_{e,k}}} = {\begin{matrix}{\lambda - \left( {\rho + {H_{k}{\hat{\Sigma}}_{k{{k - 1}}}H_{k}^{T}}} \right)} & {{- H_{k}}{\hat{\Sigma}}_{k{{k - 1}}}H_{k}^{T}} \\{{- H_{k}}{\hat{\Sigma}}_{k{{k - 1}}}H_{k}^{T}} & {\lambda - \left( {{- {\rho\gamma}_{f}^{2}} + {H_{k}{\hat{\Sigma}}_{k{{k - 1}}}H_{k}^{T}}} \right)}\end{matrix}}} \\{= {{\lambda^{2} - {\left( {{2H_{k}{\hat{\Sigma}}_{k{{k - 1}}}H_{k}^{T}} + {\rho\;\varrho}} \right)\lambda} - {\rho^{2}\gamma_{f}^{2}} + {\rho\;\varrho\; H_{k}{\hat{\Sigma}}_{k{{k - 1}}}H_{k}^{T}}} = 0}}\end{matrix} & (69) \\{\lambda_{i} = \frac{\Phi \pm \sqrt{\Phi^{2} - {4\rho\;\varrho\; H_{k}{\hat{\Sigma}}_{k{{k - 1}}}H_{k}^{T}} + {4\rho^{2}\gamma_{f}^{2}}}}{2}} & (70)\end{matrix}$where,Φ=2H _(k){circumflex over (Σ)}_(k|k−1) H _(k) ^(T) +ρl,l=1−γ_(f) ²If the following expression (71) is satisfied, one of the twoeigenvalues of the matrix R_(e,k) is positive and the other is negative,and the matrixes R_(k) and R_(e,k) have the same inertia. Therefore, theexistence condition of equation (22), i.e. (71) is obtained by using thefollowing equation (72). Here, the calculation of H_(k)K{tilde over ()}_(k) requires the same number of multiplications as O(N).−4ρlH _(k){circumflex over (Σ)}_(k|k−1)H_(k) ^(T)+4ρ²γ_(f) ²>0  (71)

$\begin{matrix}{{H_{k}{\hat{\Sigma}}_{k{{k - 1}}}H_{k}^{T}} = \frac{H_{k}{\overset{˘}{K}}_{k}}{1 - {H_{k}{\overset{˘}{K}}_{k}}}} & (72)\end{matrix}$

INDUSTRIAL APPLICABILITY

According to the present invention, as described above, the fastreal-time identification of time-invariant and time-variant systems canbe implemented by using the fast algorithm (fast H_(∞) filteringalgorithm) for the modified H_(∞) filters developed based on the newH_(∞) evaluation criterion. In addition, according to the presentinvention, the present algorithm includes, as a particular case, thefast Kalman filtering algorithm, and a term corresponding to thecovariance of system noise which is dominant in the tracking performanceof a time-varying system can be theoretically determined. Further,according to the present invention, a fast time-varying systemidentification method can be provided, which is very effectiveparticularly when a system (impulse response) is varied discontinuouslywith time, such as an echo canceller for a time-varying system whichvaries extremely as sudden line switching. Furthermore, according to thepresent invention, a system identification method can be provided, whichis applicable to echo cancellers in communication systems and acousticsystems, sound-field reproduction, and noise control.

Although embodiments of the invention have been shown and described, itis to be understood that various modifications and substitutions, aswell as rearrangements of method steps and equipment, can be made bythose skilled in the art without departing from the novel spirit andscope of the invention.

1. A system identification apparatus for performing a fast real-timeidentification of a time-invariant or time-variant system, comprising afilter robust against disturbances, said filter is formed by setting, asan H_(∞) evaluation criterion, a maximum energy gain from disturbancesweighted by ρ=1−γ_(f) ⁻² and Σ_(wi)=γ_(f) ⁻²P^_(i+1|i) to a filter errorto be smaller than a predetermined upper limit γ_(f) ², wherein saidfilter, for a state-space model as in the following equations (7) to(9), satisfies the evaluation criterion by expression (10) and is givenby the following equations (11) to (15),{circumflex over (x)} _(k+1|k+1) ={circumflex over (x)} _(k|k) +K_(s,k+1)(y _(k+1) −H _(k+1) {circumflex over (x)} _(k|k)) where,{circumflex over (x)}_(k|k): The estimate of state x_(k) at time k,obtained by using observation signals y₀ to y_(k) y_(k): Observationsignal K_(s,k+1): Filter gain H_(k): Observation matrixx _(k+1) =x _(k) +w _(k), w_(k), x_(k)εR^(N)  (7)y _(k) =H _(k) x _(k) +v _(k), y_(k), v_(k)εR  (8)z_(k)=H_(k)x_(k), z_(k)εR, H_(k)εR^(1×N)  (9) $\begin{matrix}{{\sup\limits_{x_{0},{\{\omega_{i}\}},{\{\upsilon_{i}\}}}\frac{\sum\limits_{i = 0}^{k}{{e_{f,i}}^{2}/\rho}}{{{x_{0} - {\overset{\bigvee}{x}}_{0|{- 1}}}}_{\sum\limits_{0}^{- 1}}^{2} + {\sum\limits_{i = 0}^{k}{\omega_{i}}_{\sum\limits_{\omega_{i}}^{- 1}}^{2}} + {\sum\limits_{i = 0}^{k}{{\upsilon_{i}}^{2}/\rho}}}} < \gamma_{f}^{2}} & (10)\end{matrix}${hacek over (z)}_(k|k)=H_(k){circumflex over (x)}_(k|k)  (11){circumflex over (x)} _(k+1|k+1) ={circumflex over (x)} _(k|k) +K_(s,k+1)(y _(k+1) −H _(k+1) {circumflex over (x)} _(k|k)) Filterequation  (12)K _(s,k+1) ={circumflex over (P)} _(k+1|k) H _(k+1) ^(T)(H _(k+1){circumflex over (P)} _(k+1|k) H _(k+1) ^(T)+ρ)⁻¹ Filter gain  (13)$\begin{matrix}{{{\hat{P}}_{k + {1{k}}} = {{\hat{P}}_{{{k}k} - 1} - {{{\hat{P}}_{k{{k - 1}}}\left\lbrack {H_{k}^{T}H_{k}^{T}} \right\rbrack}{R_{c,k}^{- 1}\begin{bmatrix}H_{k} \\H_{k}\end{bmatrix}}{\hat{P}}_{k{{k - 1}}}} + \sum\limits_{\omega_{k}}}}\;} & {{Riccati}\mspace{14mu}{equation}\mspace{14mu}(14)}\end{matrix}$ where, $\begin{matrix}{{e_{j,i} = {{\hat{z}}_{i{i}} - {H_{i}x_{i}}}}{R_{e,k} = {R_{k} + {\begin{bmatrix}H_{k} \\H_{k}\end{bmatrix}{{\hat{P}}_{k{{k - 1}}}\left\lbrack {H_{k}^{T}H_{k}^{T}} \right\rbrack}}}}{{R_{k} = \begin{bmatrix}\rho & 0 \\0 & {- {\rho\gamma}_{f}^{2}}\end{bmatrix}},{\Sigma_{wk} = {\gamma_{f}^{- 2}{\hat{P}}_{k + {1{k}}}}}}{{{{\hat{P}}_{k{{k - 1}}}^{- 1} + {H_{k}^{T}H_{k}}} > 0},{{\hat{P}}_{1{0}} = {ɛ_{0}I}},{ɛ_{0} > 0}}{{{0 < \rho} = {{1 - \gamma_{f}^{- 2}} \leq 1}},{\gamma_{f} > 1}}} & (15)\end{matrix}$ and the existence condition is given by the inequalityP^_(k|k−1) ⁻¹+H_(k) ^(T)H_(k)>0 which appears in (15), where, thenotation is used as follows, x_(k): State vector or just state, unknownand to be estimated, x₀: Initial state, unknown, w_(k): System noise,unknown, v_(k): Observation noise, unknown, y_(k): Observation signal,known and input to a filter, z_(k): Output signal, unknown, H_(k):Observation matrix, known, x^_(k|k): State value of the state x_(k) attime k, estimated by using observation signals y₀ to y_(k), given by thefilter equation, x^_(0|0): Initial estimate of a state, essentiallyunknown but set to 0 for convenience, K_(s,k+1): Filter gain, obtainedby matrix P^_(k+1|k,) Σ_(wk): Corresponds to the covariance matrix ofthe system noise, known in theory but unknown in advance, P^_(k|k−1):Corresponds to the covariance matrix of the error of x^_(k|k−1), givenby a Riccati equation, P^_(1|0): Corresponds to the covariance matrix ofan error in the initial state, essentially unknown but set to ε₀I forconvenience.
 2. A system identification apparatus according to claim 1,setting initial conditions of a recursive equation of the gain matrixK_(k), the auxiliary variables A_(k), S_(k), D_(k), and the stateestimate x^_(k|k), respectively, as follows,${K_{0} = 0},{A_{- 1} = 0},{S_{- 1} = \frac{\rho}{ɛ_{0}}},{D_{- 1} = 0},{{\hat{x}}_{0|0} = 0}$recursively determining the auxiliary variables A_(k) and S_(k) at timek as follows, $\begin{matrix}{{\overset{\sim}{e}}_{k} = {c_{k} + {C_{k}A_{k - 1}}}} \\{A_{k} = {A_{k - 1} - {K_{k}W_{k}{\overset{\sim}{e}}_{k}}}} \\{e_{k} = {c_{k} + {C_{k}A_{k}}}} \\{S_{k} = {{\rho\; S_{k - 1}} + {e_{k}^{T}W_{k}{\overset{\sim}{e}}_{k}}}}\end{matrix}\mspace{79mu}\begin{matrix}{\varepsilon\; R^{2 \times 1}} \\{\varepsilon\; R^{N \times 1}} \\{\varepsilon\; R^{2 \times 1}} \\{\varepsilon\; R}\end{matrix}$ obtaining a second gain matrix in which a row includingthe auxiliary variables is added to the gain matrix K_(k) as follows,${\overset{\Cup}{K}}_{k} = {\begin{bmatrix}{S_{k}^{- 1}e_{k}^{T}} \\{K_{k} + {A_{k}S_{k}^{- 1}e_{k}^{T}}}\end{bmatrix} \in R^{{({N + 1})} \times 2}}$ dividing the second gainmatrix K^(U) _(k) and obtaining first and second divisional gainmatrixes as follows, ${{\overset{\bigvee}{K}}_{k} = {{\begin{bmatrix}m_{k} \\\mu_{k}\end{bmatrix}\mspace{14mu} m_{k}} \in R^{N \times 2}}},{\mu_{k} \in R^{1 \times 2}}$determine D_(k), and obtaining a gain matrix K_(k+1) at time k+1, andobtaining a filter gain K_(s,k+1) at time k+1 by the following equationsincluding the first and second divisional gain matrixes,η_(k) =c _(k−N) +C _(k+1) D _(k−1)D _(k) =[D _(k−1) −m _(k) W _(k)η_(k)][1−μ_(k) W _(k)η_(k)]⁻¹K _(k+1) =m _(k) −D _(k)μ_(k){tilde over (K)} _(k+1)(i)=ρK _(k+1)(i, 1), i=1, . . . , NK _(s,k+1) =G _(k+1) ⁻¹ {tilde over (K)} _(k+1) , G _(k+1)=ρ+γ_(f) ⁻² H_(k+1) {tilde over (K)} _(k+1) in which η_(k)εR^(2×1), D_(k)εR^(N×1),K_(k+1)εR^(N×2), K_(s,k+1)εR^(N×1), 0<ρ=1−γ_(f) ⁻²≦1, γ_(f)>1, updatingthe filter equation (12) according to the obtained filter gain K_(s,k+1)as follows,{circumflex over (x)} _(k+1|k+1) ={circumflex over (x)} _(k|k) +K_(s,k+1)(y _(k+1) −H _(k+1) {circumflex over (x)} _(k|k)) and repeatingeach of the above means with the time being put forward, wherein theexistence of the filter is checked by using the following equation as anexistence condition suitable for fast processing, with the computationalcomplexity of O(N),−l{circumflex over (Ξ)}i+ργ _(f) ²>0, i=0, . . . , k where,${\varrho = {1 - \gamma_{f}^{2}}},{{\hat{\Xi}}_{i} = \frac{H_{i}{\overset{\sim}{K}}_{i}}{1 - {H_{i}{\overset{\sim}{K}}_{i}}}}$.
 3. A system identification apparatus according to claim 2, wherein anecho canceller is implemented by applying the filter to obtain the stateestimate x^_(k|k), producing a quasi echo as in the following equation,and canceling an actual echo by the obtained quasi echo,${\hat{d}}_{k} = {{H_{k}{\hat{x}}_{k|k}} = {\sum\limits_{i = 0}^{N - 1}{{{\hat{h}}_{i}\lbrack k\rbrack}u_{k - i}}}}$H _(k)=[u_(k) , . . . , u _(k−N+1]) where, {circumflex over (d)}_(k)Quasi-echo u_(k) Received signal N Tap number ĥ_(i)[k] Estimate ofimpulse response of echo path.
 4. A system identification apparatusaccording to claim 1, wherein an echo canceller is implemented byapplying the filter to obtain the state estimate x^_(k|k), producing aquasi echo as in the following equation, and canceling an actual echo bythe obtained quasi echo,${\hat{d}}_{k} = {{H_{k}{\hat{x}}_{k|k}} = {\sum\limits_{i = 0}^{N - 1}{{{\hat{h}}_{i}\lbrack k\rbrack}u_{k - i}}}}$H _(k) =[u _(k) , . . . , u _(k−N+1]) where, {circumflex over (d)}_(k)Quasi-echo u_(k) Received signal N Tap number ĥ_(i)[k] Estimate ofimpulse response of echo path.